#ifndef __MATRIX_OPS_H__
#define __MATRIX_OPS_H__
#include <iostream>

namespace algebra
{
	template<typename T> T det(const Matrix<T,3,3> & m)
	{
		T result = m.at(0,0)*m.at(1,1)*m.at(2,2) + m.at(1,0)*m.at(2,1)*m.at(0,2) + m.at(0,1)*m.at(1,2)*m.at(2,0);
		result -= m.at(2,0)*m.at(1,1)*m.at(0,2) + m.at(1,0)*m.at(0,1)*m.at(2,2) + m.at(2,1)*m.at(1,2)*m.at(0,0);
		return result;
	}

	template<typename T> T det(const Matrix<T,2,2> & m)
	{
		return m.at(0,0)*m.at(1,1) - m.at(1,0)*m.at(0,1);
	}

	template<typename T1, typename T2, int sizex1, int sizey1, int sizex2,int sizey2>
		bool areOfSameDimension(const Matrix<T1,sizex1,sizey1> &m1, const Matrix<T2,sizex2,sizey2> m2)
		{
			return ( (sizex1 == sizex2) && (sizey1 == sizey2));
		}

	template<typename T> 
		Matrix<T,2,2> inverse(const Matrix<T,2,2> & m)
		{
			T detvalue = det(m);
			Matrix<T,2,2> result;
			result.at(0,0) = m.at(1,1);
			result.at(1,0) = -m.at(1,0);
			result.at(0,1) = -m.at(0,1);
			result.at(1,1) = m.at(0,0);
			return result/detvalue;
		}

	template<typename T, int square_side> 
		bool isInvertible(const Matrix<T,square_side,square_side> & m)
		{
			return det(m) != 0;
		}

	template<typename T>
		Matrix<T,1,3> operator^(const Matrix<T,1,3> & v1, const Matrix<T,1,3> &v2)
		{
			Matrix<T,1,3> result;
			result.at(0,0) = v1.at(0,1) * v2.at(0,2) - v1.at(0,2) * v2.at(0,1);
			result.at(0,1) = v1.at(0,2) * v2.at(0,0) - v1.at(0,0) * v2.at(0,2);
			result.at(0,2) = v1.at(0,0) * v2.at(0,1) - v1.at(0,1) * v2.at(0,0);
			return result;
		}
	
	template<typename T,int side>
		Matrix<T,1,side> normalize(const Matrix<T,1,side> & v1)
		{
			T sum = 0;
			for(int i=0;i<side;++i)
				sum+=v1[i] * v1[i];
			sum = (T)(1.0/sqrt(sum));
			Matrix<T,1,side> result;
			for(int i=0;i<side;++i)
				result[i] = v1[i] * sum;
			return result;
		}
	
	template<typename T, int square_side> 
		Matrix<T,square_side,square_side> transpose(const Matrix<T,square_side,square_side> & m)
		{
			Matrix<T,square_side,square_side> result;
			for(int i=0;i<square_side;++i)
				for(int j=0;j<square_side;++j)
				{
					result.at(i,j) = m.at(j,i);
				}
			return result;
		}
}

#endif

